Ejercicio: Datos satelitales en GRASS GIS
Trabajamos con imágenes Sentinel 2 en GRASS GIS
Datos Sentinel 2
- Lanzamiento: Sentinel-2A en 2015, Sentinel-2B en 2017
- Tiempo de revisita: ~5 días
- Cobertura sistemática de áreas terrestres y costeras entre los 84°N y 56°S
- 13 bandas espectrales con resolución espacial de 10 m (VIS y NIR), 20 m (red-edge y SWIR) y 60 m (otras)
ESA - Satélites Copernicus Sentinel. Más información en: https://www.copernicus.eu/en/about-copernicus/infrastructure/discover-our-satellites
Distribución de bandas de Sentinel 2 comparadas con Landsat
Extensiones para datos Sentinel
- i.sentinel.download: descarga productos Copernicus Sentinel de Copernicus Open Access Hub
- i.sentinel.import: importa datos Sentinel descargados de Copernicus Open Access Hub
- i.sentinel.preproc: importa y realiza corrección atmosférica y topográfica de imágenes S2
- i.sentinel.mask: crea máscaras de nubes y sombras para imágenes S2
Ver Sentinel wiki para más detalles)
Recientemente, se sumaron nuevos miembros en la familia i.sentinel:
i.sentinel.coverage: comprueba la cobertura de área de las escenas de S1 o S2 seleccionadas
i.sentinel.parallel.download: descarga imagenes Sentinel en paralelo
Para conectarse al Copernicus Open Access Hub a través de i.sentinel.download, se necesita ser usuario registrado
Crear el archivo
SENTINEL_SETTING.txten el directorio$HOME/gisdata/con el siguiente contenido:
your_username
your_password
Niveles de procesamiento Sentinel 2
- L1C: Reflectancia a tope de atmósfera o Top of Atmosphere (TOA). Disponibles desde el lanzamiento.
- L2A: Reflectancia Superficial o Bottom of Atmosphere (BOA), i.e., los datos han sido corregidos para remover los efectos de la atmósfera. Sólo desde 2019.
Archivo de datos Sentinel
Long Term Archive (LTA)
Todos los productos (1C o 2A) de más de un año son movidos fuera de línea y se requiere un tiempo de espera para ponerlos a disposición del usuario. Esto dificulta la automatización de tareas con productos de más de 12 meses de antigüedad.
Iniciar GRASS GIS, crear nuevo mapset y establecer región computacional
import subprocess
import sys
# Ask GRASS GIS where its Python packages are to be able to start it from the notebook
sys.path.append(
subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)
# Importar los paquetes python de GRASS
import grass.script as gs
import grass.jupyter as gj
# Iniciar GRASS
session = gj.init(grassdata, location, mapset)Crear un nuevo mapset sentinel2
Definir la región computacionalal radio urbano de Córdoba
Búsqueda y descarga de datos S2
Instalar la extensión i.sentinel
Lista de escenas disponibles que intersectan la región computacional
Lista de escenas disponibles que contienen la región computacional
Descargar la escena seleccionada - NO EJECUTAR
Como la descarga desde el Copernicus Open Access Hub toma su tiempo, vamos a descargar la escena Sentinel 2 que usaremos y moverla a HOME/gisdata/s2_data
Hagamos una prueba con datos del LTA…
Importar datos Sentinel 2 a GRASS GIS
1. Importar con corrección atmosférica: i.sentinel.preproc
Productos nivel 1C
Para obtener un valor de AOD, tenemos 2 opciones:
A. Estimar el valor desde un grafico
B. Descargar un archivo y el valor sera estimado
Obtener AOD de
http://aeronet.gsfc.nasa.gov
- Estación ARM_Cordoba o Pilar_Cordoba
- Seleccionar fechas de inicio y final
- Seleccionar:
Combined fileyAll points - Descargar y descomprimir (el archivo final tiene extensión .dubovik)
- Pasar el archivo con la opción
aeronet_file
Mapa de elevación
- r.in.srtm.region: importa (y re-proyecta) los mosaicos SRTM que cubren la región computacional, parchea los mosaicos e interpola datos faltantes
- r.in.nasadem: importa (y re-proyecta) los mosaicos de NASADEM que cubren la región computacional y parchea los mosaicos
Si el DEM es más chico que la región computacional, sólo la región cubierta por el DEM será corregida atmosféricamente…
Ejemplo
2. Importar sin corrección atmosférica (as is): i.sentinel.import
Productos nivel 2A
Imprimir información sobre las bandas antes de importarlas
Importar bandas seleccionadas, recortar y reproyectar al vuelo
Listar bandas importadas y revisar metadatos
Balance de colores y composiciones
Asignar grey como paleta de colores
Ajuste de colores para una composición RGB color natural
Mostrar la combinación RGB 432
Tarea
Realizar balance de colores y mostrar combinacion falso color NIR-RED-GREEN
Máscara de nubes y sombras de nubes
Identificar y enmascarar nubes y sus sombras
# identify and mask clouds and clouds shadows: i.sentinel.mask
i.sentinel.mask -s --o \
blue=T20JLL_20200330T141049_B02_10m \
green=T20JLL_20200330T141049_B03_10m \
red=T20JLL_20200330T141049_B04_10m \
nir=T20JLL_20200330T141049_B08_10m \
nir8a=T20JLL_20200330T141049_B8A_20m \
swir11=T20JLL_20200330T141049_B11_20m \
swir12=T20JLL_20200330T141049_B12_20m \
cloud_mask=cloud \
shadow_mask=shadow \
scale_fac=10000 \
mtd=$HOME/gisdata/s2_data/S2B_MSIL2A_20200330T141049_N0214_R110_T20JLL_20200330T182252.SAFE/GRANULE/L2A_T20JLL_A016009_20200330T141532/MTD_TL.xmlVisualización de la salida: Nubes y sombras de nubes identificadas con i.sentinel.mask
Índices de agua y vegetación
Definir región computacional
Establecer máscara
Estimación de los índices de vegetación
Instalar extensión i.wi
Estimación de índice de agua
Visualización de los resultados
Segmentación
Instalar la extensión i.superpixels.slic
Listar los mapas y crear grupos y subgrupos
Ejecutar i.superpixels.slic
Convertir el resultado a vector
Ejecutar i.segment
Convertir el resultado a vector
Mostrar NDVI junto con las 2 salidas de la segmentación
Tarea
Ejecutar cualquiera de los 2 métodos de segmentación con diferentes parámetros y comparar los resultados
Clasificación supervisada
Tarea
- digitalizar áreas de entrenamiento para 3 clases con g.gui.iclass
- guardarlas en un mapa vectorial:
training
Clasificación supervisada con Maximum Likelihood
Convertir el vector de áreas de entrenamiento a raster
Generar archivos de firma espectral
Realizar la clasificación por Maximum Likelihood
Añadir etiquetas a las clases
Clasificación supervisada con Maximum Likelihood
Clasificación supervisada con Machine Learning
Instalar la extensión r.learn.ml
Realizar la clasificación por RF
Añadir etiquetas a las clases
Clasificación supervisada con Random Forest
Tarea
Comparar los resultados de ambos tipos de clasificación supervisada a través del índice Kappa
Hay un módulo r.kappa
Post-procesamiento y validación
- r.reclass.area para eliminar pequeñas áreas, enmascarar nuevos valores y rellenar los huecos con r.neighbors o r.fillnulls
- convertir la salida en vector y ejecutar v.clean con
tool=rmarea - r.kappa para la validación (idealmente también digitalizar una muestra de prueba)